A van der Waals density functional mapping of attraction in DNA dimers 
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The dispersion interaction between a pair of parallel DNA double-helix structures is investigated 
by means of the van der Waals density functional (vdW-DF) method. Each double-helix structure 
consists of an infinite repetition of one B-DNA coil with 10 base pairs. This parameter-free density 
functional theory (DFT) study illustrates the initial step in a proposed vdW-DF computational 
strategy for large biomolecular problems. The strategy is to first perform a survey of interaction 
geometries, based on the evaluation of the van der Waals (vdW) attraction, and then limit the eval- 
uation of the remaining DFT parts (specifically the expensive study of the kinetic-energy repulsion) 
to the thus identified interesting geometries. Possibilities for accelerating this second step is detailed 
in a separate study. For the B-DNA dimer, the variation in van der Waals attraction is explored at 
relatively short distances (although beyond the region of density overlap) for a 360° rotation. This 
study highlights the role of the structural motifs, like the grooves, in enhancing or reducing the 
vdW interaction strength. We find that to a first approximation, it is possible to compare the DNA 
double strand at large wall-to-wall separations to the cylindrical shape of a carbon nanotube (which 
is almost isotropic under rotation). We compare our first-principles results with the atom-based 
dispersive interaction predicted by DFT-D2 [J. Comp. Chem. 27, 1787 (2006)] and find agreement 
in the asymptotic region. However, we also find that the differences in the enhancement that occur 
at shorter distances reveal characteristic features that result from the fact that the vdW-DF method 
is an electron-based (as opposed to atom-based) description. 



I. INTRODUCTION 

The key characteristics of a working biomolecular sys- 
tem is an enormous richness of structural complexity 
and a robust working principle, molecular recognition, 
for identifying geometries that optimize the intermolec- 
ular binding and alignment. Stronger covalent or ionic 
binding determines the structure inside the molecules 
and provides resilience towards statistical fluctuations.^ 
This permanence is, for example, of utmost importance in 
the preservation of the information con tained in our de- 
oxyribonucleic acid (DNA) genome.l^ The binding that 
is of relevance for life processes, e.g., the molecular- 
recognition matching^HSl gf genes, is weaker to allow the 
reversible operations that Nature needs.'^' The molecu- 
lar recognition arises as a delicate balance between steric 
hindrance, electrostatics, and van der Waals (vdW) at- 
traction, the latter also termed the London dispersion 
interaction. 

The search for a deeper understanding of such 
biomolecular operation motivates development of a 
parameter-free computational theory that both pro- 
vides transferability and computational efficiency.^ ^ The 
structural complexity of the biomolecular systems im- 
plies that the supramolecular system may express itself 
in a multitude of ways, and it is not certain that a given 
empirical (or semi-empirical) interaction model remains 
applicable for all emerging configurations and for varying 
charging states. 

Density Functional Theory (DFT) is a highly valued 
condensed-matter theory tool that is our work horse in 
predictive computational theory of traditional materials 
problems.'^ In such systems, the electron density remains 
high between the atoms, for example, in a bulk structure. 



DFT also works excellently for individual molecules but 
it has until recently lacked an account of the truly nonlo- 
cal correlations that underpin vdW interactions between 
constituents in a molecular system. The issue is that 
molecular systems are inherently sparseJ^^HlIl they must 
also contain low electron-density regions between the 
molecular fragments. The sparseness is even the defining 
quality when Nature puts molecular recognition to work 
among biomolecules. Nevertheless, the last decade have 
seen developments that position DFT to overcome this 
previous limitation. 

The vdW density functional (vdW-DF) methocPE^ 
has a track recorcpSl as a general-purpose nonempirical 
DFT that can characterize interactions in sparse and 
soft biomolecular systems. The vdW-DF method ac- 
counts for vdW interactions by introducing true nonlo- 
cality in the density functional. Being built as a physics- 
based description and within a constraint-based design, 
it has the electron-based description that helps trans- 
ferability. This is, for example, important for systems 
where it is essential to describe the vdW binding or at- 
traction for several typical separations.E^HlIl This is true 
even if the plasmon-type description may not capture 
the full complexity in the description of the far-apart- 
regime for some systems that effectively behave as a 
low-dimensional metal.E^lEil There exist efficient imple- 
mentations also for self-consistent evaluations. One such 
algorithnp2uses a fast Fourier transform approach which 
accelerates evaluation for medium-to-large size systems 
while real-space evaluation approaches'^ 23 26 j^^^y hold 
advantages for very-large-scaled systems, as discussed in 
Ref.Ull 

This paper reports that a nonempirical vdW-DF char- 
acterization of the variation in the vdW attraction among 
biopolymers is indeed feasible. In vdW-DF, this vdW 



attraction is reflected in a nonlocal correlation term 
E^^. Other terms, including a single-particle kinetic en- 
ergy term T^ and the semilocal parts of the exchange- 
correlation energy i?°^ also contribute to the vdW-DF 
total energy E'^^^~^^ but (for objects that are effectively 
neutral) it is the vdW attraction, contained in i?"\ that 
dominates when the density overlap between fragments 
can be ignored. Our vdW-DF characterization presented 
here thus computes the £^"' variation to map the mor- 
phology dependence in the vdW attraction between a 
pair of parallel, infinitely-repeated double-stranded (ds) 
DNA helices. We find that the vdW interaction of such 
DNA dimer is sensitive to the alignment of the ma- 
jor DNA structural motifs, particularly different orienta- 
tions of the major and minor grooves of the two ds-DNA 
macromolecules. We also identify and discuss a set of sys- 
tematic changes that arise in the vdW attraction when 
the problem is investigated across a range of interaction 
distances.'i^^^ We further compare the i<^"' variation to 
the evaluation of the semiempirical vdW term of a vdW- 
extended DFT method,'^!! namely DFT-D2.'22l We use and 
extend analytical interaction resultP^lEll that apply for 
the far (but not asymptotically far) regime, to illustrate 
the difference in nature between our electron-based vdW- 
DF study and the atom-based DFT-D approach. More- 
over, we contrast this interpretation of the vdW attrac- 
tion in a DNA dimer with corresponding results for a 
dimer of carbon nanotubes (CNTs). 

With this study we propose a vdW-DF computational 
strategy for large-biomolecules, the strategy begins non- 
empirical studies of the interactions by a mapping of the 
vdW-attraction term £■"'. This is suggested to optimize 
computational efficiency. We observe that the i?"' map- 
ping of vdW attraction can today be carried out at the 
cost of about ten wall time minutes per DNA dimer in- 
teraction geometry. That is much smaller than the time 
it takes to evaluate the variation in the kinetic energy 
Ts because of an excellent scaling of real-space evalua- 
tions of i?"', Ref. [TTI The cost of computing Tg through 
wavefunctions is the same for DFT-D and vdW-DF so 
our suggested strategy could also b e rele vant, for exam- 
ple, for large-system DFT-D studies .f^SESxhe motivation 
for our proposal is that a complete nonempirical vdW- 
DF evaluation can be focused on the relevant interaction 
morphologies that the E'^^ mapping identifies. The possi- 
bilities for using an adapted Harris schemeSSEHl to accel- 
erate the complete vdW-DF evaluation is discussed and 
assessed in a closely related pubhcation, Ref. 18 

The paper is organized as follows. The ds-DNA and 
DNA dimer systems are described and discussed in Sec- 
tion II, which also introduces the particular CNT struc- 
ture investigated here. In Section III we present the com- 
putational method used. Section IV consists of a study of 
the nonlocal correlation energy variation as the distance 
between the DNAs is changed and the two structures are 
rotated. It also contains comparisons of our vdW-DF- 
based results in the far (nearly asymptotic) region to the 
dispersion interaction predicted by DFT-D2,'23 and with 




FIG. 1: Along-axis and side views of a minimal repeat-unit 
cell model of the ds-DNA macromolecule. The atoms are P 
(yellow), O (dark red), H (cyan), N (blue), and C (gray). 
Figure created using XCrySDen.-^ 



results of similar studies for CNTs. 



II. GENOME STRUCTURE, A DNA MODEL 

The primary structure of DNA is made up by four dif- 
ferent nucleobases, the adenine (A), guanine (G), cyto- 
sine (C) and thymine (T) bases, each covalently linked to 
a sugar. These building blocks are joined to a strand by 
phosphate groups connecting the sugars. The two DNA 
strands form a double helix with (predominantly) hydro- 
gen bonds linking base pairs, thus forming the DNA sec- 
ondary structure. The repeated pairing between a purine 
base (A or G) and a smaller pyrimidine (T or C) produces 
a constant diameter for the double helix. 

One of us was previously involved in a vdW-DF study 
of the nucleobases, focusing on the stacking interactions.'^ 
Even if the vdW forces by themselves may have limited 
selectivity (for small molecules like base pairs), the vdW 
forces are essential in the building of this secondary struc- 
ture because they drive, optimize, and hold a structural 
assembly until more stable bonds can be formed. This 
mechanism is fundamental in the case of molecular recog- 
nition and in DNA replication. The previous study pro- 
vides vdW-DF analysis and results for the Rise and Twist 
of the base pairs inside (that is, between the backbones 
of) the ds-DNA structure. The nonlocal correlation con- 
tribution E^^ to vdW-DF will by itself have difficulties in 
aligning small molecules but the vdW forces still position 
these components (base pairs) close enough that the cor- 
rugation from kinetic-energy repulsion is expressed. The 
previous vdW-DF study^ reported a very good agreement 
between vdW-DF characterizations and experimental in- 
formation (for example, correlations in base-pair rotation 
angles as extracted from the Nucleic Acid Database^). 
This progress highlights possibilities of a broader appli- 
cation of the vdW-DF method to soft matter challenges 
and in particular to refine the description of the vdW 
attraction between biopolymers. 

We focus our large-scale vdW-DF interaction study 
on dimers of a ds-DNA model structure in the B-form. 



Both ds-DNA fragments are assumed to have the period- 
ically repeated sequence of nucleobases GCAATACGGT. 
The ds-DNA atomic structure is built from idealized 
coordinates P^H^ In our model the axes of the two ds-B- 
DNA fragments are straight and parallel to each other. 

Figure [T] shows the atomic positions in the minimal- 
repeat unit cell for one ds-B-DNA molecule. On the most 
coarse level it reflects a cylinder-like form for DNA but 
it also shows both major and minor grooves of the dou- 
ble helix. This DNA model contains 635 atoms and is 
limited to one coil, but it is periodic and is used as the 
repeat unit cell of a DFT calculation so that the DNA is 
still described as a macromolecular system (through the 
infinite repetition of the unit cell along the DNA axis). 

The polar sugar-phosphate strand on the external part 
of the chain can form favorable interactions with ions 
in a solvent. Each phosphate group then carries a neg- 
ative charge which is balanced by positive ions, the 
counter ions, in the surrounding liquid. The counter 
ions stabilize the structure, making the DNA-ion system 
charge neutral. Correlation between the concentrations 
in the counter ions may themselves contribute further 
to an attraction between the DNA-ion system and an- 
other organic molecule. However, we shall not consider 
such effects here. Instead we investigate the vdW bind- 
ing arising from nonlocal electron correlation effects be- 
tween two ds-DNA structures assuming (and enforcing in 
our DFT calculations) that the DNA repeated-unit cell 
model structure itself remains charge neutral. A forth- 
coming study addresses the effects on the vdW attraction 
of the charging by counter ions. 

The double-helix structure is very regular, and when 
viewed along the axis appears almost isotropic (left panel 
of Fig. [I]) . This average isotropy of DNA simplifies our 
study of the sensitivity of the vdW binding to macro- 
molecular structural motifs. The average isotropy also 
helps in the comparison with CNT dimer interactions 
as we can define an approximate "radius" of DNA. We 
can thus compare the vdW interactions as a function of 
wall-to-wall distances A in dimers of DNA and dimers 
of CNT, using the analysis of Refs. [27l80H 32l For DNA 
we use the radius that is defined as the average of the 
helix backbone P and O atoms distances from the DNA 
center (8.5 A and 9.7 A). This yields the approximate 
DNA radius (r)DNA ~ 9.1 A (shown by a double-arrow 
line in Figure IT]), close to that of the standard literature 
description where the diameter is taken to be ~ 20 A. For 
a given center-to-center separation d of the DNA dimer 
structure we thus use A = d— 2(r)DNA as the wall-to-wall 
separation of the DNA dimer. 

The approximate similarity of DNA with the cylinder 
form of a CNT represents one major structural motif of 
the structure. Another major structural motif is obvi- 
ously the existence of grooves (right panel), a feature 
which is distinctly different from what characterizes the 
CNTs. 

The systems studied here are sparse^^ so that the vdW 
interaction (with electrostatics) dominates in the inter- 




FIG. 2: Left: Example of a system obtained starting from a 
single 33.8 A period of the ds-DNA that is copied, translated 
and rotated inside an enlarged unit cell. Right: Portion of 
the final system studied. Here we have piled a number of 
unit cells (not all shown) on top of each other in preparation 
of the real-space evaluation of Eq. S. Figure created using 
XCrySDen.I^ 



molecular regions of low electron density. By basing the 
function on weaker forces, Nature ensures a truly remark- 
able resilience of our genetic code.l'Sl 



III. COMPUTATIONAL DETAILS 

DFT calculations of interacting systems are effectively 
limited by the computational challenge of accurately cal- 
culating the kinetic energy, a step which requires solving 
the noninteracting particle Schrodinger equation in three 
dimensions. One could fear that the nonlocal nature of 
the correlations that define the dispersive forces, i.e., the 
vdW-DF evaluation of the i?"', may also represent a com- 
putational challenge for large systems. However, efficient 
algorithms are now in place and there are, effectively, no 
other computational bottlenecks for biomolecular vdW- 
DF and/or DFT-D studies, like those presented here, 
than the evaluation of the (noninteracting) kinetic en- 
ergy Ts. 

We note that since i?"' is far less costly (for large sys- 
tems) than any calculation of eigenstates, it is wise to 
begin a first-principle biomolecular interaction study by 
mapping out the _E"' variation; this saves the costly de- 
termination of the kinetic-energy repulsion to relevant 
interaction morphologies. Here we pursue this first map- 
ping step of evaluating E^^ to identify optimal DNA- 
dimer configurations. We also characterize the sensitivity 
of the vdW binding to variations in alignment of DNA 
structural motifs. 

Figure ^ illustrates our process in our .©"'-mapping. 
It is formally the first step in an overall Harris-type ap- 
proach for an accelerated vdW-DF evaluation of inter- 
molecular interactions.'iS 

The DNA-dimer electron density n^; is needed for eval- 
uating i?"' [ud] for the dimer. For the individual DNA we 



calculate the electron density in a Linear Combination of 
Atomic Orbitals (LCAO) characterization, available in 
the DFT code GPAWpi' As an approximation of n^ we 
use the superposition 
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of two copies of the individual density, n'^^ . The super- 
script indicates a standard self-consistent GGA calcula- 
tion. 

The calculation of the electron density that makes up 
71^2^^ is carried out at the F point. The length of a period 
along the ds-DNA is 33.8 A, and for determining the 
electron density we use an orthorhombic unit cell of size 
39.5 X 40.4 X 33.8 A^. The electron density is described 
on a grid with a spacing of approximately 0.20 A. 

The extended system for the superposition n^ is cre- 
ated by first enlarging the unit cell of the individual 
DNA by factors 2.3 x 2.3 x 1 to obtain a unit cell of size 
90.9 X 92.9 X 33.8 A^ (in Figure [5] the original unit cell 
is enlarged only along the first axis, for visual reasons). 
Then the density for the second ds-DNA is added, appro- 
priately translated and/or rotated around the DNA axis. 
For rotated ds-DNA the values of the electron density 
on the spatial grid are interpolated from the grid values 
before rotation. 

The set of leftmost panels of Fig. [2] summarizes this 
step in our vdW-DF study. By also rotating the ds-DNA 
densities we can map out various mutual alignments of 
the groove structure and study the ramifications for the 
vdW attractions (as described in vdW-DF or in DFT-D). 

For each of the considered separations d of the dimer 
and for each of the relative angles of the fragments we 
evaluate ii^"' by carrying out the integral 



Ec\n] = - / n(ri)0[n](ri,r2)n(r2)dridr2 (2) 



1 
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for the superposition density n^ of Eq. (IT]). Here (j)[n\ 
is given from a tabulated kernel function, but (/)[n] still 
refiects the overall density variation in the interacting 
components.'i^ElMl Our code for the postprocess calcu- 
lation of the E^^ value is a real-space code. For the pe- 
riodicity along the DNA axes we perform the real-space 
evaluation step by replicating the DNA density along the 
z direction (the DNA axis) , for a total length of up to 540 
A (depending on the need in each calculation) and then 
evaluate ^ with one spatial coordinate restricted to the 
central unit cell, as described in Ref. [17] and illustrated 
also in Ref. |331 This is indicated by the right panel of 
Fig. [2] 

The ds-DNA structure can be roughly approximated 
by a CNT-like cylindrical structure (although filled) 
when viewed at large wall-to-wall separation. For a com- 
parison we present CNT calculations using a (15,15) CNT 
with the electron density (of the individual CNT) ob- 
tained in a 40.7 x 40.7 x 2.46 A^ unit ceU that contains 60 
atoms. The radius of the (15,15) CNT is (r)cNT=10.0 A, 
approximately the size of the DNA radius. For the CNT 



we have chosen to use the GPAW finite differences mode 
(instead of LCAO) as a basis for the individual CNT 
electron density variation. We also set a grid for the 
description of the electron density that has a spacing 
less than 0.135 A between grid points, and we use a 
Monkhorst-Pack fc-point sampling of the Brillouin zone 
with a 1x1x16 mesh. Like for the DNA dimer, the CNT 
dimer density is obtained as a superposition of the indi- 
vidual densities. 

In an additional set of comparisons, we also report 
DFT-D2 calculations of the vdW attraction for dimers 
of DNA and for dimers of CNTs. These are based on the 
dimer atomic configurations. The DFT-D calculations 
are carried out as follows. For each of the DNA/CNT 
atoms in the original unit cell, we sum the pair contri- 
butions over atoms in the repeated copies of the DNA or 
CNT unit cell, with the repetitions extending over 473 
A (or 8890 atoms) for the DNA-dimer problem and over 
404 A (or 9900 atoms) for the CNT-dimer problem. 



IV. RESULTS AND DISCUSSION 

A. A mapping of mutual vdW attraction in a soft 
twinning of two B-DNA coils 

We first explore the vdW sensitivity of the DNA inter- 
actions as one DNA is rotated (rigidly) around another 
DNA, as illustrated in the insert of Fig.[3J In the rotation 
6 the individual DNA is not rotated around its own axis. 
Instead the situation corresponds to keeping the DNA 
dimer axes fixed in space and rotating the left hand side 
DNA by the angle —9 around its own axis while simulta- 
neously rotating the right hand side DNA around its own 
axis by the angle —9. Effectively we are thus exploring 
what variation arises in the vdW attraction effects when 
twinning the two B-DNA coils. We restrict this twinning 
to situations that create no significant density overlap. 

Each of the curves in Fig. |3] corresponds to a certain 
relative twinning rotation 9 and shows the interaction 
E^^ at various dimer separations (wall-to-wall distance A 
or center-to-center distance d) from values smaller than 
expected at the dimer binding, A = 2.94 A, and up to 
A = 12.4 A. 

We note that the value of _E"' varies strongly with wall- 
to-wall separation A, but except at the very small sepa- 
ration A = 2.94 A, the order of the curves of largest and 
smallest contribution remain the same. 

We cover the variation in the rotation angle 9 up to 
90° with steps of 10°, and we find variation in the DNA 
vdW binding 9. At about A — 3.7 A (within binding- 
distance range) the variation, per unit length, in E'"' with 
9 amounts to about 7 meV/A. 

There are several local minima in the variation when 
viewed at fixed distance > 3 A (at 6^ = 0°, 40°, and 70°). 
The variation in the attraction with 9 is correlated with 
the alignment of the grooves, but the presence of several 
local minima also shows that detailed structure in the 
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FIG. 3: Nonlocal correlation energy per length _B" jL of the 
DNA dimer as a function of rotation angle Q and separation. 
Both the center-to-center d and wall-to-wall distances A are 
shown in the plot on the bottom and top horizontal axis, 
respectively. The picture in the insert shows the geometry of 
the system as the angle Q is varied. 



DNA further influences the interaction. For example, the 
DNA is not a continuous double helix, but consists of a 
series of base pairs and the associated parts of the back 
bone, at 3.4 A apart along the DNA axis, and the choice 
of bases in the base pairs varies along the DNA. Thus 
we should not expect smooth, largely monotonous curves 
when the relative orientations of the DNA fragments are 
changed. 



B. Interaction effects by the alignment of motifs 

The curves in Figure [3] show significant variation in 
£""' with relative orientation of the DNA fragments. We 
therefore further explore the sensitivity to alignment of 
structural motifs by keeping one of the two DNA frag- 
ments fixed and rotating the other by an angle 0. Our 
study creates a 360° mapping of Ef- jL versus 0, the rota- 
tion of one B-DNA coil. The results are shown in Figure 
fThe interaction is evaluated at a separation A = 4.5 
which is slightly larger than the expected binding dis- 
tance. At closer distances some orientations will result 
in overlap of high densities at the outer O atoms on the 
two fragments (because DNA is not truly cylindrical). 

Fig. HI demonstrates a significant and rapid variation of 
the vdW binding with the alignment of motifs, i.e., the 
alignment or disalignment of the major and minor groves. 
There is an approximate symmetry of the nonlocal cor- 
relation energy around the central drop (at 183°), with 
the value —14.2 meV/A, or 12.0 meV/A lower in energy 
than the least attractive orientations at that separation, 
namely at 125° and 235°. 

The Ef^ (or vdW attraction) minimum at 183° corre- 
sponds to the alignment of the two structures in a way 
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FIG. 4: The vdW binding as a function of the rotation of 
one of the DNA molecules around its axis while keeping the 
other at a fixed orientation. The separation of the dimer is 
A = 4.5 A. 



that, for each of the fragments considered, two of the 
outermost O atoms face each other directly. This also 
means that the minor and major groove of the first struc- 
ture are aligned with the corresponding grooves on the 
second DNA, as depicted in the inset picture associated 
to this point. 

The two orientations (at = 125° and 235°) with the 
largest E^ jL values have a vdW attraction of only —2.4 
meV/A, respectively —2.6 meV/A, compared to the sit- 
uation of the two DNA fragments being far apart. This 
orientation has a large amount of vacuum between the 
two structures along their entire length. 

Two other (minimum) features correspond to a clear 
enhancement in attraction, even if more modest. These 
local-minimum angles appear as approximately symmet- 
ric with respect to the global minimum in the ^ variation. 
They appear at 40° and 328° with vdW attraction -6.7 
meV/A and with —7.9 meV/A. This is 4-6 meV/A lower 
in energy compared to the least attractive orientation. 
These two local minima correspond to an alignment of 
the two DNA in which only one of the external O atoms 
faces another O atom in the other DNA. 

Altogether, our vdW-DF evaluation thus makes it pos- 
sible to detect variations in the nonlocal correlation en- 
ergy within a ~12 meV/A broad range due to the align- 
ment of the structural motifs. 



C. A comparison with vdW-extended DFT and 
with a CNT dimer 

We continue the analysis by comparing and contrast- 
ing our vdW-DF results with results from a summation 
of pair contributions, using a traditional atom-based dis- 
persive interaction form, here as provided by Grimme 
for the dispersion term D of DFT-D in the version DFT- 
D2.29 Figure[5]repeats the 9 = 40° and 9 = 50° curves for 
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FIG. 5: DNA and CNT vdW interaction per length as a 
function of wall-to-wall distance A. 



the DNA vdW interaction and includes also the curves 
for the dispersion term D of DFT-D for the same two 
rotations. 

For separations larger than about A = 6 A, we find 
reasonable agreement between the results of the vdW-DF 
and DFT-D procedures (Fig. [s]) . At smaller separations 
we expect the results to differ because at close range it 
becomes important that the dispersion interaction arises 
mainly in the valence electron region, not at the atomic 
centers, as assumed in DFT-D. Indeed, we do see a differ- 
ence in the results at small separations, with diminished 
attraction in the results from DFT-D compared to those 
from the vdW-DF method. 

We also compare in Figure [5] to the attraction of a pair 
of parallel (15,15) CNTs that have approximately the 
same radius as DNA, (r)cNT = 10 A. The CNT data are 
calculated with vdW-DF and DFT-D. The attraction of 
the CNTs is stronger than for the DNA dimer. We note 
that there are about 30% more atoms per length (rele- 
vant for DFT-D, although also the species of the atoms 
matter) and more integrated electron density per length 
(relevant for vdW-DF) i n the (15,15) CNT than in the 
DNA. In a previous studjEMl gf dimers of polyethylene 
(PE), polypropylene (PP), and polyvinyl chloride (PVC) 
two of us found a strong dependence of the vdW interac- 
tion (in a simplified method) on the integrated amount 
of electron charge (per length) in the polymer. In a full 
vdW-DF calculation, including _E"', of PE in dimers and 
in a crystal, we founcP^ the vdW term _E"' to depend 
strongly on the separation of the molecules and less, but 
not insignificantly, on the relative orientation. 

We stress that there are significantly larger differences 
between the estimates of the vdW attraction for a CNT 
dimer based on our vdW-DF evaluation and the DFT- 
D2 evaluation, than for a DNA dimer in Fig. [6j Some 
of the difference must be ascribed to the difference that 
exists in the underlying GPAW calculation of the single- 
fragment electron density (a courser-grained LCAO pro- 
cedure for the DNA and a finite-different approach for 
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FIG. 6: The vdW interaction in DNA and CNT dimers using 
vdW-DF and DFT-D. 



CNT). However, we expect (and are presently testing in a 
continuation vdW-DF study) that the differences mostly 
reflect the fact that for vdW-DF it is important not only 
how many electrons are in a molecule but also how they 
are distributed in space.B2l It is also the case that for 
vdW-extended DFT methods, (which use parameters to 
specify the per-atom contribution to the vdW interaction 
strength) there is now an awareness that sp2-hybridized 
carbon (the form that roughly applies for the CNTs) will 
need a different paramctrization than carbon in other or- 
ganic molecules (like DNA). 



D. On the validity of a near-asymptotic interaction 

form 



Fig. [6] analyzes the asymptotic forms of the vdW-DF 
and DFT-D2 results for the (15,15)-CNT and DNA in- 
teraction curves for the 9 = 40° orientation. The largest 
center-to-center distance d = A -1- 2(r) is here approxi- 
mately 60 A, i.e., at the furthest separation dmax the ra- 
tio to the DNA radius is dmax/('')DNA ~ 6. This means 
that only the points furthest out can possibly be con- 
sidered truly asymptotic. In the asymptotic and almost- 
asymptotic range the = 40° and 50° curves are identi- 



cal and we therefore only show the — 40° curves in this 
plot. 

Beyond d w 50 A and up to dmax, which is a small 
range here, the interactions for the DNA and the CNT, 
both from vdW-DF and DFT-D calculations, scale ap- 
proximately as — d~P where p = 5. This is the asymp- 
totic scaling expected for parallel cylinders. The slope of 
—d~^ is drawn in the top panel of Figure [6| as a guiding 
line. In general the interaction of parallel cylinders at 
equal radius h at large (but not necessarily asymptotic) 
separation d is given by the generalized hypergeometric 
function j,F2, as described in the Appendix. This applies 
for cylinders of a continuous material, which our sys- 
tems with atoms and varying electron distribution arc 
not quite. 

In the bottom panel of Figure |6] we plot the interac- 
tion per length, times d^ . In the truly asymptotic region 
for a pair of cylinders we expect a vanishing slope, and 
for separations a bit smaller we expect the next terms 
in the expansions ( |A.4 ) for filled cylinders and (A.6) for 
hollow cylinders to contribute. As seen in the bottom 
panel of Figure [6j neither the DNA nor the CNT dimers 
are in the fully asymptotic region at d < dmax (the 
curves are not totally flat). However, after adding the 
first few terms of the expansions (A.4) or (A.6) we find 



good agreement with the form of the generalized hyper- 
geometric function 3i^2 outside the near-binding region, 
A = d-2(r) < 4-10 A. 

Whether a pair of cylinders are filled or hollow does not 
affect the lowest order term in the asymptotic expansion, 
the (i~^-term. However, it does change by a factor of 2 
the coefficient on the next order term, the d~^-term, and 
thus the interaction curve in the not-quite-asymptotic 
region. When the DNAs or CNTs are close to each other 
the separation is too small for the expressions in terms of 
the generalized hypergeometrical functions to be valid. 

It is clear that DNA is not a regular cylinder structure 
and that the existence of grooves must play an essential 
role, as also seen for the interaction as a function of rel- 
ative orientation of the DNA dimer. However, far from 
the binding region, basing the description on an assump- 
tion of the electrons being distributed in a filled cylinder 
is a good approximation. 

The analytical evaluation refiects the morphology of 
the interaction fragments and is therefore a good approx- 
imation to the pair-potential summation that underpins 
a DFT-D evaluation. Fig.[6]shows that the analytical de- 
scription applies well as an approximation to the DFT- 
D results when the fragments are sufficiently removed 
to also make the above-stated assumptions meaningful. 
However, Fig. p^ also shows that our E'^'^-hdsed calcula- 
tion of the vdW attraction maintains differences for the 
near-asymptotic form out to further distances than does 
the DFT-D2-based evaluation. 

We ascribe these differences to the vdW-DF ability to 
include multipole and some collectivity effects through its 
plasmon-pole descriptioiP^l an d its emphasis on reflecting 
the electron-density variation.E^EZl 



E. Towards a full vdW-DF interaction study 

We note that at distances A < 4 A it is also impor- 
tant to include the remaining components of the vdW-DF 
method, see Ref. [ini 

In a closely related study^^ we assess the possibility 
of using an adaption of the Harris schem#2t26l to accel- 
erate the evaluation also of the remaining vdW-DF (or 
DFT-D) components, predominantly seeking to bypass 
repeated evaluations of the kinetic energy variation T^. 
This study shows that a high degree of accuracy can be 
achieved. A forthcoming paper will report the results for 
the DNA dimer problem and will detail effects that the 
DNA charging state might have on the vdW attraction. 



V. SUMMARY AND OUTLOOK 

In this paper, we show that the vdW-DF functional 
has the potential to describe extended systems with the 
accuracy of DFT, thus opening a way to the description 
of complex soft phenomena. 

We find that our computational strategy for first- 
principle vdW-DF characterization (with our moderate- 
level access to HPC) can give some qualitative and even 
semi-quantitative results for biomolecular systems even 
before the full vdW-DF (or for example DFT-D) evalua- 
tion of the kinetic-energy repulsion effects is completed. 

An accelerated vdW-DF computational frameworWi^ 
study can be expected to work well for investigations of 
vdW bound systems. Here, we have illustrated the pos- 
sibility of an even faster initial mapping by evaluating 
the vdW-binding component. An additional component 
in the strategy is to adapt the ideas of the Harris scheme 
as explored in a parallel study^ The motivation such 
strategy is easily stated: even if our DNA model struc- 
ture is only a model system of the full secondary genome 
structure, the individual ds-DNA still contains 635 atoms 
and the dimer system contains 1270 atoms. The size of 
the problem prevents us (at the supercomputing facilities 
to which we have access) from performing an ordinary 
potential-energy calculation by means of a conventional 
DFT implementation. More realistic studies will contain 
even more atoms per unit cell, or have other similar com- 
plications. Thus any computational simplification able 
to keep most of the original accuracy is welcome. That 
observation is true whether or not we wish to pursue a 
vdW-DF or a DFT-D study 

The nonlocal correlation energy _E"', which contains 
the essential information about the vdW binding, is ac- 
cessible with a limited computational effort. We have 
found that this S"' evaluation has essentially perfect 
scaling up to at least 2000 cores in our real space 
implementation.121 

We are thus proposing to initialize vdW-DF studies 
of large biomolecular interaction problems by first map- 
ping out what interaction geometries are plausible, from 
knowing the variation in the vdW attraction. By begin- 



ning the vdW-DF identification of optimal interaction 
geometries with this vdW-attraction step we are pursu- 
ing a course that is similar to that Nature uses in its own 
molecular-recognition search. 



Here i? is a (positive) prefactor that reflects the suscep- 
tibility of the material in the cylinders and whose value 
can be set by investigating the asymptotic form. 

We rewrite this expression with the generalized hyper- 
geometric function a^Pj of one variable 
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Appendix: Scaling of interaction at large separations 

With some robust assumptions it is possible to pro- 
vide analytical results for the vdW interaction per length, 
Ea/L, for two parallel (infinitely repeated) cylindrical 
structures when these are far, but not asymptotically 
far, apart .122121132] rpj^^ assumptions are that the matter 
is continuous, a slightly simplified plasmon-pole approx- 
imation, and that the (screened) effective susceptibility 
can be approximated as having equivalent magnitudes 
for response along the cylinder axis and tangent. There 
are some differences in this analytical evaluation for hol- 
low (relevant for carbon nanotubes) and filled (relevant 
for DNA) cylinders but in both cases the results can be 
expressed in terms of a generalized hypergeometric func- 
tion. For either types of interacting systems, we consider 
two cylinders of equal radius b at center-to-center dis- 
tance d. 

For a pair of massive cylinders a two- variable expres- 
sion is given in the book by Mahanty and NinhanpSBZI jj^ 
terms of Appells hypergeometric function F4 
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which is the expansion used in the bottom panel of Figure 
[6] for the pair of DNA (approximated as filled cylinders) . 

For a pair of infinitely thin, hollow cylinders, where the 
electron density can be described by a radial 5-function, 
the interaction i s give n by a similar generalized hyperge- 
ometric functiorP322l 



T^hollow 



64 



15 5-1 1-4^ 
2'2'2' ' ' d2 



(A.5) 



The result is also valid for a slightly more general 
choice of susceptibility tensors. '^^ ^^ For the function 



(A.5 1 the expansion in b/d yields 
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This expansion is used in the bottom panel of Figure|6]for 
the pair of CNT (approximated as pipes, that is, hollow 
cylinders). 



* Electronic address: |schroder@chalmers.se| Corresponding 
author 



^ E. Schrodinger, What is life? (Cambridge University 
Press, Cambridge, 1967). 



9 



^ J.D. Watson and F.H.C. Crick, Nature 171, 4356 (1953). 
^ H.M. Berman, W.K. Olson, D.L. Beveridge, J. West- 
brook, A. Gelbin, T. Demeny, S.-H. Hsieh, A.R. Srini- 

vasan, and B. Schneider, Biophys. J. 63, 751 (1992); 

[http://ndbserver.rutg ers.edu| 
■* M.C. Williams, L. J. Maher III, eds., "Biophysics of DNA- 

Protein Interactions" (Springer, New York, 2011). 
^ A.L. Lehninger, D.L. Nelson, and M.M. Cox, Principles of 

Biochemistry, (Worth Publishers, Inc., New York, 1993). 
^ For example, S.H. Gellman, Chem. Rev. 97, 1231 (1997) 

and associated collected reviews in Chem. Rev. 97; P.J. 

Edmonson et al, Int. J. Mod. Sci. 9, 154 (2008). 
''' V.R. Cooper, T. Thonhauser, A. Puzder, E. Schroder, B.I. 

Lundqvist, and D.C. Langreth, J. Amer. Chem. Soc. 130, 

1304 (2008). 
* S. Li, V.R. Cooper, T. Thonhauser, B.I. Lundqvist, and 

D.C. Langreth, J. Phys. Chem. B 113, 11166 (2009). 
^ K. Burke, J. Chem. Phys. 136, 150901 (2012). 
^° H. Rydberg, N. Jacobson, P. Hyldgaard, S.I. Simak, 

B.I. Lundqvist, and D.C. Langreth, Surf. Sci. 532—535, 

606 (2003). 
" D.C. Langreth, M. Dion, H. Rydberg, E. Schroder, P. 

Hyldgaard, and B.I. Lundqvist, Intern. J. Quant. Chem. 

101, 599 (2005). 
^^ D.C. Langreth, B.I. Lundqvist, S.D. Chakarova-Kack, 

V.R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, 

J. Kleis, L. Kong, S. Li, P.G. Moses, E. Murray, A. Puzder, 

H. Rydberg, E. Schroder, and T. Thonhauser, J. Phys.: 

Cond. Matter 21, 084203 (2009). 
^^ M. Dion, H. Rydberg, E. Schroder, D.C. Langreth, and 

B.I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004); 95, 

109902(E) (2005). 
^^ T. Thonhauser, V.R. Cooper, S. Li, A. Puzder, P. 

Hyldgaard, and D.C. Langreth, Phys. Rev. B 76, 125112 

(2007). 
^^ K. Lee, E.D. Murray, L. Kong, B.I. Lundqvist, and 

D.C. Langreth, Phys. Rev. B 82, 081101(R) (2010). 
1® K. Berland and P. Hyldgaard, J. Chem. Phys. 132, 134705 

(2010). 
^^ K. Berland, 0. Borck, and P. Hyldgaard, Comp. Phys. 

Commun. 182, 1800 (2011). 
^* K. Berland, E. Londero, E. Schroder, and P. Hyldgaard, 

A Harris-type van der Waals density functional scheme, 

arXiv: 1303.3762 [cond-mat.mtrl-sci]. 
^^ Y.U. Barash, Fiz. Tverd. Tela. 30, 1578 (1988); B. E. Ser- 

nelius and P. Bjork, Phys. Rev. B 57, 6592 (1998). 
^° J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 

073201 (2012); S. Lebegue, J. Harl, T. Gould, J. Angyan, 

G. Kresse, and J.F. Dobson, Phys. Rev. Lett. 105, 196401 

(2010). 
^^ J. P. Perdew, J. Tao, P. Hao, A. Ruzsinszky, G.I. Csonka, 

and J.M. Pitarke, J. Phys.: Condens. Matter 24, 424207 

(2012); A. Ruzsinszky, J. P. Perdew, J. Tao, G.I. Csonka, 

and J. M. Pitarke, Phys. Rev. Lett. 109, 233203 (2012). 
^^ G. Roman-Perez and J.M. Soler, Phys. Rev. Lett. 103, 

096102 (2009). 
^^ E. Ziambaras, J. Kleis, E. Schroder, and P. Hyldgaard, 

Phys. Rev. B 76, 155425 (2007). 



^'' A. Gulans, M.J. Puska, and R.M. Nieminen, Phys. Rev. B 
79, 201105(R) (2009). 

^^ D. Nabok, P. Puschnig, C. Ambrosch-Draxl, Comp. Phys. 
Commun. 182, 1657 (2011). 

^^ Open-source tool JuNoLo for real-space and fast-fourier 
transform non-selfconsistent vdW-DF total-energy eval- 
uation; P. Lazic, N. Atodiresei, M. Alaei, V. Caciuc, 
S. Bliigel, and R. Brako, Comp. Phys. Comm. 181, 371 
(2010). 

^■^ J. Kleis, E. Schroder, and P. Hyldgaard, Phys. Rev. B 77, 
205422 (2008). 

^* X. Wu, M.C. Vargas, S. Nayak, V. Lotrich, and 
G. Scoles, J. Chem. Phys. 115, 8748 (2001); S. Grimme, 
J. Comp. Chem. 25, 1463 (2004). 

2^ S. Grimme, J. Comp. Chem. 27, 1787 (2006). 

3° E. Schroder and P. Hyldgaard, Surf. Sci. 532, 880 (2003). 

^^ E. Schroder and P. Hyldgaard, Mater. Sci. Engin. C 23, 
721 (2003). 

^^ J. Kleis, P. Hyldgaard, and E. Schroder, Comp. Mater. Sci. 
33, 192 (2005). 

^^ J. Harris, Phys. Rev. B 31, 1770 (1985). 

^^ W.M.C. Foulkes and R. Haydock, Phys. Rev. B 39, 12520 
(1989). 

35 V.K. Nikuhn, Zh. Tekhn. Fiz. XLI, 41 (1971) [Sov. Phys. 
- Techn. Phys. 16, 28 (1971)]. 

'^^ R.G. Gordon and Y.S. Kim, J. Chem. Phys. 56, 3122 
(1972). 

^"^ A. Kokalj Comp Mater. Sci. 28, 155 (2003). Code avail- 
able from |http://www.xcryscien. org/ 1 

3* S. Arnott and D.W.L. Hukins, Biochem. Biophys. Research 
Commun. 47, 1504 (1972). 

3^ S. Arnott, D.W.L. Hukins, and S.D. Dover, Biochem. Bio- 
phys. Research Commun. 48, 1392 (1972). 

"^^ The idealized atomic coordinates for the here-investigated 
base-pair sequence in our infinitely repeated ds-DNA 
model system were put together by Prof. Jason Kahn, 
University of Maryland. The coordinates are available 
as part of his instructive rasmol tutorial, available at 
jhttp://www. biochem. umd.edu/biochem/kahn/teachjes] 

^^ Open-source, grid-based PAW-method DFT code GPAW, 
http://wiki.fysik.dtu.dk/gpaw/ 

*^ K. Berland and P. Hyldgaard, An analysis of van der 
Waals density functional components: Binding and corru- 
gation of benzene and C60 on boron nitride and graphene, 
arXiv: 1303.0389 [cond-mat.mes-hall]. 

■^^ K. Berland, S.D. Chakarova-Kack, V.R. Cooper, D.C. Lan- 
greth, and E. Schroder, J. Phys.: Cond. Matt. 23, 135001 
(2011). 

■^'^ J. Kleis and E. Schroder, J. Chem. Phys. 122, 164902 
(2005). 

'^'' J. Kleis, B.I. Lundqvist, D.C. Langreth, and E. Schroder, 
Phys. Rev. B 76, 100201(R) (2007). 

'^^ J. Mahanty and B.W. Ninham, Dispersion Forces (Aca- 
demic Press, London, 1976). 

''^ Please note that there is a typographical error in Ref. |46] 
in what corresponds to the last expression below. 



